library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
#library(corrplot)
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('keep.trailing.zeros',TRUE)
source("C:/Users/jtame/Documents/GitHub/RISKPLOTS/CODE/auxplots.R")
RRPLOTS and flchain
odata <- flchain
odata$chapter <- NULL
pander::pander(table(odata$death))
rownames(odata) <- c(1:nrow(odata))
data <- as.data.frame(model.matrix(Surv(futime,death)~.,odata))
data$`(Intercept)` <- NULL
dataFL <- as.data.frame(cbind(time=odata[rownames(data),"futime"],
status=odata[rownames(data),"death"],
data))
pander::pander(table(dataFL$status))
dataFL$time <- dataFL$time/365
Exploring Raw Features with RRPlot
convar <- colnames(dataFL)[lapply(apply(dataFL,2,unique),length) > 10]
convar <- convar[convar != "time"]
topvar <- univariate_BinEnsemble(dataFL[,c("status",convar)],"status")
pander::pander(topvar)
topv <- min(5,length(topvar))
topFive <- names(topvar)[1:topv]
topFeature <- RRPlot(cbind(dataFL$status,dataFL[,topFive[1]]),
title=topFive[1])


par(op)
## With Survival Analysis
RRanalysis <- list();
idx <- 1
for (topf in topFive)
{
RRanalysis[[idx]] <- RRPlot(cbind(dataFL$status,dataFL[,topf]),
timetoEvent=dataFL$time,
atRate=c(0.90,0.80),
title=topf)
idx <- idx + 1
par(op)
}












names(RRanalysis) <- topFive
Reporting the Metrics
pander::pander(t(RRanalysis[[1]]$keyPoints),caption="Threshold values")
Threshold values
| Thr |
73.000 |
69.000 |
68.000 |
53.000 |
50.00000 |
| RR |
4.013 |
4.399 |
4.537 |
5.592 |
1.00000 |
| RR_LCI |
3.740 |
4.045 |
4.160 |
4.406 |
0.00000 |
| RR_UCI |
4.305 |
4.783 |
4.947 |
7.096 |
0.00000 |
| SEN |
0.581 |
0.713 |
0.733 |
0.967 |
1.00000 |
| SPE |
0.883 |
0.790 |
0.776 |
0.216 |
0.00329 |
| BACC |
0.732 |
0.752 |
0.755 |
0.591 |
0.50164 |
ROCAUC <- NULL
CstatCI <- NULL
RRatios <- NULL
LogRangp <- NULL
Sensitivity <- NULL
Specificity <- NULL
for (topf in topFive)
{
CstatCI <- rbind(CstatCI,RRanalysis[[topf]]$c.index$cstatCI)
RRatios <- rbind(RRatios,RRanalysis[[topf]]$RR_atP)
LogRangp <- rbind(LogRangp,RRanalysis[[topf]]$surdif$pvalue)
Sensitivity <- rbind(Sensitivity,RRanalysis[[topf]]$ROCAnalysis$sensitivity)
Specificity <- rbind(Specificity,RRanalysis[[topf]]$ROCAnalysis$specificity)
ROCAUC <- rbind(ROCAUC,RRanalysis[[topf]]$ROCAnalysis$aucs)
}
rownames(CstatCI) <- topFive
rownames(LogRangp) <- topFive
rownames(Sensitivity) <- topFive
rownames(Specificity) <- topFive
rownames(ROCAUC) <- topFive
pander::pander(ROCAUC)
| age |
0.822 |
0.810 |
0.833 |
| kappa |
0.682 |
0.667 |
0.696 |
| lambda |
0.665 |
0.650 |
0.680 |
| creatinine |
0.589 |
0.574 |
0.605 |
pander::pander(CstatCI)
| age |
0.775 |
0.774 |
0.764 |
0.785 |
| kappa |
0.671 |
0.671 |
0.659 |
0.684 |
| lambda |
0.657 |
0.657 |
0.644 |
0.671 |
| creatinine |
0.585 |
0.585 |
0.572 |
0.599 |
pander::pander(LogRangp)
| age |
0.00e+00 |
| kappa |
4.90e-175 |
| lambda |
4.41e-145 |
| creatinine |
2.67e-67 |
pander::pander(Sensitivity)
| age |
0.581 |
0.558 |
0.602 |
| kappa |
0.319 |
0.298 |
0.340 |
| lambda |
0.300 |
0.279 |
0.321 |
| creatinine |
0.288 |
0.269 |
0.309 |
pander::pander(Specificity)
| age |
0.883 |
0.873 |
0.892 |
| kappa |
0.900 |
0.891 |
0.908 |
| lambda |
0.899 |
0.890 |
0.907 |
| creatinine |
0.865 |
0.854 |
0.875 |
meanMatrix <- cbind(ROCAUC[,1],CstatCI[,1],Sensitivity[,1],Specificity[,1])
colnames(meanMatrix) <- c("ROCAUC","C-Stat","Sen","Spe")
pander::pander(meanMatrix)
| age |
0.822 |
0.775 |
0.581 |
0.883 |
| kappa |
0.682 |
0.671 |
0.319 |
0.900 |
| lambda |
0.665 |
0.657 |
0.300 |
0.899 |
| creatinine |
0.589 |
0.585 |
0.288 |
0.865 |
Train Test Set
trainsamples <- sample(nrow(dataFL),0.5*nrow(dataFL))
dataFLTrain <- dataFL[trainsamples,]
dataFLTest <- dataFL[-trainsamples,]
pander::pander(table(dataFLTrain$status))
pander::pander(table(dataFLTest$status))
Cox Modeling
ml <- BSWiMS.model(Surv(time,status)~.,data=dataFLTrain,loops=0)
sm <- summary(ml)
pander::pander(sm$coefficients)
Table continues below
| age |
0.0985 |
0.0727 |
0.0797 |
0.086 |
0.718 |
0.608 |
| flc.grp |
0.1042 |
0.0635 |
0.0908 |
0.115 |
0.605 |
0.714 |
| creatinine |
0.4552 |
0.3074 |
0.4103 |
0.590 |
0.548 |
0.719 |
Table continues below
| age |
0.725 |
0.741 |
0.630 |
0.75 |
0.19084 |
0.86869 |
| flc.grp |
0.725 |
0.631 |
0.740 |
0.75 |
0.01121 |
0.31941 |
| creatinine |
0.725 |
0.549 |
0.745 |
0.75 |
0.00344 |
0.00935 |
| age |
27.02 |
26.162 |
-0.12016 |
1 |
| flc.grp |
5.51 |
8.543 |
-0.01039 |
1 |
| creatinine |
3.10 |
0.245 |
-0.00525 |
1 |
Test results
index <- predict(ml,dataFLTest)
rtestdata <- cbind(dataFLTest$status,ppoisGzero(index,h0))
rrAnalysisTest <- RRPlot(rtestdata,atRate=c(0.90,0.80),
timetoEvent=dataFLTest$time,
title="Test: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






By Risk Categories
obsexp <- rrAnalysisTest$OERatio$atThrEstimates
expObs(obsexp,"Test: Expected vs. Observed")

pander::pander(obsexp)
| Total |
988 |
927 |
1052 |
1333 |
0.741 |
0.696 |
0.789 |
4.85e-23 |
| low |
260 |
229 |
294 |
322 |
0.808 |
0.713 |
0.913 |
4.40e-04 |
| 90% |
157 |
133 |
184 |
167 |
0.943 |
0.801 |
1.102 |
4.86e-01 |
| 80% |
574 |
528 |
623 |
845 |
0.680 |
0.625 |
0.738 |
6.35e-23 |
Time to Event Analysis
rrAnalysisdata <- rrAnalysisTest
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction:
rrAnalysisdata$timetoEventData$eTime and
rrAnalysisdata$timetoEventData$expectedTime
| 566654 |
0 * * * |
two.sided |
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction:
rrAnalysisdata$timetoEventData$eTime[highrisk] and
rrAnalysisdata$timetoEventData$expectedTime[highrisk]
| 265800 |
8.99e-57 * * * |
two.sided |
timesdata <- expObsTime(rrAnalysisdata,title="Test: Expected vs Observed")

pander::pander(timesdata)
Table continues below
| 1Q |
10.5 |
6.99 |
| Median |
12.4 |
11.30 |
| 3Q |
13.2 |
12.91 |
Table continues below
| 1Q |
2.93 |
26.5 |
| Median |
6.96 |
44.0 |
| 3Q |
11.02 |
70.7 |
| 1Q |
9.81 |
2.35 |
| Median |
11.31 |
4.05 |
| 3Q |
12.65 |
6.11 |
Cox Calibration
op <- par(no.readonly = TRUE)
calprob <- CoxRiskCalibration(ml,dataFLTrain,"status","time")
( 15.13201 , 17550.98 , 1.272783 , 971 , 1068.674 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
The RRplot() of the calibrated model
index <- predict(ml,dataFLTrain)
calrdata <- cbind(dataFLTrain$status,ppoisGzero(index,calprob$h0))
rrAnalysisCalTrain <- RRPlot(calrdata,atRate=c(0.90,0.80),
timetoEvent=dataFLTrain$time,
title="Cal. Train: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=calprob$timeInterval)






By Risk Categories
obsexp <- rrAnalysisCalTrain$OERatio$atThrEstimates
expObs(obsexp,"Cal: Expected vs. Observed")

pander::pander(obsexp)
| Total |
971 |
911 |
1034 |
842 |
1.15 |
1.08 |
1.23 |
1.29e-05 |
| low |
270 |
239 |
304 |
220 |
1.23 |
1.09 |
1.38 |
9.40e-04 |
| 90% |
139 |
117 |
164 |
116 |
1.20 |
1.01 |
1.42 |
3.63e-02 |
| 80% |
562 |
516 |
610 |
505 |
1.11 |
1.02 |
1.21 |
1.27e-02 |
Time to Event Analysis
rrAnalysisdata <- rrAnalysisCalTrain
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction:
rrAnalysisdata$timetoEventData$eTime and
rrAnalysisdata$timetoEventData$expectedTime
| 248652 |
0 * * * |
two.sided |
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction:
rrAnalysisdata$timetoEventData$eTime[highrisk] and
rrAnalysisdata$timetoEventData$expectedTime[highrisk]
| 193150 |
2.01e-08 * * * |
two.sided |
timesdata <- expObsTime(rrAnalysisdata,title="Cal: Expected vs Observed")

pander::pander(timesdata)
Table continues below
| 1Q |
10.4 |
7.85 |
| Median |
12.4 |
11.85 |
| 3Q |
13.2 |
13.11 |
Table continues below
| 1Q |
3.14 |
37.2 |
| Median |
7.29 |
64.7 |
| 3Q |
11.34 |
108.6 |
| 1Q |
14.3 |
3.78 |
| Median |
16.2 |
6.30 |
| 3Q |
18.2 |
9.13 |
Checking the test set
index <- predict(ml,dataFLTest)
rtestdata <- cbind(dataFLTest$status,ppoisGzero(index,calprob$h0))
rrAnalysisCalTest <- RRPlot(rtestdata,atRate=c(0.90,0.80),
timetoEvent=dataFLTest$time,
title="Cal. Test: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=calprob$timeInterval)






By Risk Categories
obsexp <- rrAnalysisCalTest$OERatio$atThrEstimates
expObs(obsexp,"Cal Test: Expected vs. Observed")

pander::pander(obsexp)
| Total |
988 |
927 |
1052 |
818 |
1.21 |
1.13 |
1.29 |
8.66e-09 |
| low |
260 |
229 |
294 |
201 |
1.29 |
1.14 |
1.46 |
5.66e-05 |
| 90% |
157 |
133 |
184 |
104 |
1.51 |
1.28 |
1.77 |
1.11e-06 |
| 80% |
574 |
528 |
623 |
513 |
1.12 |
1.03 |
1.21 |
8.61e-03 |
Time to Event Analysis
rrAnalysisdata <- rrAnalysisCalTest
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))
Wilcoxon signed rank test with continuity correction:
rrAnalysisdata$timetoEventData$eTime and
rrAnalysisdata$timetoEventData$expectedTime
| 202579 |
0 * * * |
two.sided |
highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))
Wilcoxon signed rank test with continuity correction:
rrAnalysisdata$timetoEventData$eTime[highrisk] and
rrAnalysisdata$timetoEventData$expectedTime[highrisk]
| 173740 |
0.0606 |
two.sided |
timesdata <- expObsTime(rrAnalysisdata,title="Cal Test: Expected vs Observed")

pander::pander(timesdata)
Table continues below
| 1Q |
10.5 |
6.99 |
| Median |
12.4 |
11.30 |
| 3Q |
13.2 |
12.91 |
Table continues below
| 1Q |
2.93 |
42.4 |
| Median |
6.96 |
70.5 |
| 3Q |
11.02 |
113.2 |
| 1Q |
15.7 |
3.76 |
| Median |
18.1 |
6.49 |
| 3Q |
20.3 |
9.79 |